% Fourier法解一维扩散方程
% 参考：Biner 相场编程
% Gitee Repo

clc
clear


L=1;
dx=0.02;
dt=0.00001;

x=(-L:dx:L)';
n=size(x,1);

D=1;

if mod(n,2) == 1
	k1 = 2*pi/(2*L)*((1:ceil(n/2))'-1);
	k2 = -flip(k1);
	k2(end,:)=[];
	kx = [k1;k2];
	clear k1
	clear k2
else
	error('i cannot handle it');
end


u0 = zeros(n,1);
u1 = zeros(n,1);
fu0 = zeros(n,1);
fu1 = zeros(n,1);

u0 = exp(-10*(x-0.1*L).^2);
fu0 = (fft(u0));

for tick = 1:10000
    % 显式法
	fu1 = fu0 - dt*D*kx.^2.*fu0;
    % 隐式法
    fu1 = fu0./(1+dt*D*kx.^2);

    fu0=fu1;

    % 作为对比的有限差分法
    u1(2:n-1) = u0(2:n-1) +dt*D/(dx^2)*(u0(1:n-2)-2*u0(2:n-1)+u0(3:n));
    u0 = u1;

	if mod(tick,100) == 0
		clf
		axis([-L  L 0 1])
        axis equal
		hold on

		plot(x,real(ifft(fu1)))
		plot(x,u1)

		drawnow
		pause(0.1)
	end
end
